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Abstract 


This  report  reviews  a  numerical  model  for  calculating  the  evolution  of  a  break¬ 
ing  wave.  The  model  is  the  combination  of  a  modified  version  of  RIPPLE  which 
was  originally  developed  at  Los  Alamos  National  Laboratory  (Kothe  et  ah,  1991) 
and  the  A:  -  e  turbulence  model.  In  the  model,  finite  difference  solutions  to  the  in¬ 
compressible  Reynolds  equations  for  the  mean  flow  field  and  the  k  —  t  equations  for 
the  turbulent  field  are  obtained  on  a  nonuniform  mesh.  The  free  surface  locations 
are  represented  by  the  volume  of  fluid  (VOF)  data  on  the  mesh.  A  two-step  projec¬ 
tion  method  is  used  for  the  mean  flow  solutions,  aided  by  the  incomplete  Cholesky 
conjugate  gradient  technique  solving  the  Poisson  equation  for  the  mean  pressure 
field.  Advections  of  momentum  in  Reynolds  equations  and  turbulent  kinetic  energy 
and  dissipation  rate  in  the  k  —  e  equations  are  estimated  by  the  combination  of  the 
upwind  method  and  the  central  difference  method.  Several  numerical  examples, 
including  the  runup  and  rundown  of  nonbreaking  and  breaking  solitary  waves,  are 
given.  Agreement  between  the  experimental  data  and  the  numerical  results  is  very 
good. 


1  Introduction 

Wave  transformation  phenomena  in  shallow  water  are  of  great  theoretical  and 
practical  importance.  Nearshore  breaking  waves  are  central  to  nearly  all  coastal 
processes,  including  coastal  currents  and  sediment  transport.  Breaking  waves  are 
also  responsible  for  the  production  of  air-bubbles  and  sea-water  droplets,  which 
are  important  in  the  consideration  of  the  pollutant  transport.  From  the  viewpoint 
of  fluid  mechanics,  breaking  wave  transformation  processes  still  possess  many  chal¬ 
lenging  but  unsolved  problems. 

Intensive  researches  have  been  performed  in  the  last  two  decades  to  study  the 
detailed  mechanisms  of  breaking  waves.  Battjes  (1988)  summarized  both  the  ex¬ 
perimental  findings  and  numerical  approaches  for  the  breaking  waves  in  the  surf 
zone.  Although  many  laboratory  experiments  have  been  conducted  to  measure  the 
velocity  field  and  the  free-surface  profile  during  the  breaking  process  (e.g.  Stive, 


1980;  Mizuguchi,  1986;  Skjelbria,  1987;  Nadaokaet  al,  1989;  Ting  and  Kirby,  1994, 
1995,  and  1996),  direct  numerical  solutions  for  breaking  waves  are  rare.  Tliis  is 
primarily  due  to  the  lack  of  understanding  of  turbulent  flows  associated  with  wave 
breaking. 

The  early  approximation  of  wave  breaking  simulation  stemmed  from  the  numer¬ 
ical  computation  of  depth  averaged  equations  (DAE),  which  include  both  shallow 
water  equations  (SWE)  and  Boussinesq  equations  (Peregrine,  1967).  The  energy 
dissipation  due  to  the  breaking  process  was  incorporated  into  the  equations  through 
certain  dissipative  terms  (Abbott  et  al.,  1978;  Svendsen  et  al.,  1978).  Because 
the  dimension  is  reduced  by  one  in  the  DAE,  the  computational  expense  is  much 
cheaper  than  that  for  full  equations  and  thus  this  approach  can  be  carried  on  rather 
large  scale  simulation  such  as  tsunami  propagation  and  runup  (Liu  et  al.,  1993). 
Even  today,  it  is  still  a  very  active  research  area  of  modeling  wave  transformation 
in  surf  zone  with  the  use  of  the  DAE  (Abbott  et  al.,  1985;  Zelt,  1991;  Karambas 
and  Koutitas,  1992;  Schaffer  et  al.,  1993;  Kobayashi  and  Karjadi,  1993;  Eldeberky 
and  Battjes,  1996).  By  taking  advantage  of  the  simple  form  of  the  DAE,  how¬ 
ever,  we  should  also  realize  its  limitations.  Since  the  DAE  approach  requires  a 
single  value  of  free  surface  displacement,  this  approach  cannot  predict  the  detailed 
configuration  of  the  free  surface  during  the  breaking.  Furthermore,  this  approach 
cannot  provide  the  information  of  the  generation  and  transport  of  vorticity  and 
turbulence.  These  details  can  only  be  recovered  if  the  full  hydrodynamic  equations 
are  solved. 

Another  important  approach  to  simulate  the  initiation  of  wave  breaking  is  to 
employ  the  potential  flow  theory.  Before  a  wave  breaks,  wave  motions  are  es¬ 
sentially  irrotational,  except  in  boundary  layer  regions.  Therefore,  the  evolution 
process  can  be  described  theoretically  by  a  potential  flow  theory.  The  potential 
flow  with  a  free  surface  can  be  numerically  computed  with  the  boundary  integral 
equation  method  (BIEM).  The  earliest  pioneer  of  using  the  BIEM  to  study  the 
wave  breaking  problem  is  Longuet-Higgins  and  Cokelet  (1976).  The  method  was 
further  improved  by  Dommermuth  et  al.  (1988)  and  Grilli  et  al.  (1994).  With  the 
use  of  the  BIEM,  the  free  surface  displacement  is  no  longer  limited  to  be  single 
valued  as  in  the  DAE.  The  detailed  shape  of  a  curled  plunging  jet  can  be  simulated 
very  well  and  the  depth-dependent  velocity  information  is  also  available.  However, 
the  BIEM  has  a  severe  limitation  that  the  computation  cannot  continue  after  the 
plunging  jet  touches  the  front  free  surface.  Therefore,  many  interesting  and  im¬ 
portant  features  of  wave  breaking  such  as  the  generation  and  transport  of  vorticity 
and  turbulence  cannot  be  studied  with  the  BIEM  method. 

Free  of  all  the  limitations  discussed  above,  the  full  Reynolds  equations  can 
describe  the  mean  (ensemble  averaged)  motion  of  any  turbulent  flow.  The  Reynolds 
equations  have  the  similar  form  as  the  Navier-Stokes  equations  (NSE)  and  the 
former  is  actually  the  special  case  of  the  latter  when  the  flow  is  laminar.  The 
numerical  solution  to  the  NSE  was  studied  intensively  in  1960’s  and  1970’s.  The 
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earliest  solver  to  the  NSE  together  with  the  Marker  and  Cell  (MAC)  method,  a 
free  surface  tracking  technique,  was  proposed  by  Harlow  and  Welch  (1965).  The 
method  has  been  continuously  improved  in  the  following  decades.  One  of  the 
outstanding  methods  is  so  called  the  projection  method  developed  by  Chorin  in 
1968  and  1969.  This  method  is  very  robust  and  will  be  adapted  in  our  model, 
serving  as  the  basic  solver  to  the  Reynolds  equations. 

Another  important  issue  concerning  breaking  wave  computations  is  the  free 
surface  tracking  technique.  Considering  the  complex  geometry  of  the  free  surface 
during  the  breaking  process,  we  immediately  rule  out  the  traditional  height  func¬ 
tion  method  that  requires  the  free  surface  displacement  to  be  single  valued.  The 
MAC  method  (Harlow  and  Welch,  1965)  does  not  have  the  limitation  of  geome¬ 
try  complexity,  but  it  requires  too  much  storage  to  save  the  marker  information. 
Because  the  marker  is  in  general  not  located  at  the  place  where  the  velocity  is  de¬ 
fined,  the  movement  of  those  markers  must  be  based  on  the  interpolated  velocity 
that  may  lead  to  rather  large  accumulated  errors.  The  surface  segment  method 
(Miyata,  1986)  or  the  surface  marker  (SM)  method  (Chen  et  ah,  1991)  significantly 
reduces  the  storage  by  tracking  the  markers  on  the  free  surface  only,  but  it  intro¬ 
duces  additional  difficulties  of  reordering  surface  marker  when  the  reconnection  of 
free  surface  occurs.  Such  difficulties  may  becomes  unsolvable  when  a  strong  break¬ 
ing  takes  place  that  generates  many  small  droplets  and  bubbles.  Attempting  to 
track  the  density  change  within  each  computational  cell  instead  of  the  free  surface 
location,  the  volume  of  fluid  (VOF)  method  (Hirt  and  Nichols,  1981)  provides  a 
robust  alternative  for  updating  the  free  surface.  With  the  considerations  of  both 
accuracy  and  efficiency,  we  decide  to  adapt  the  VOF  method  in  our  current  model. 

Combining  the  solver  to  the  NSE,  such  as  the  projection  method,  and  free  sur¬ 
face  tracking  technique,  such  cis  the  MAC  or  the  VOF  method,  we  can  in  principle 
solve  any  laminar  flow  with  free  surface,  i.e.,  non  breaking  wave  (Chan  and  Street, 
1970).  For  breaking  wave,  because  of  the  generation  of  strong  turbulence,  addi¬ 
tional  turbulence  model  must  be  added  to  represent  the  effect  of  turbulence  on 
the  mean  flow  motion.  In  practice,  however,  the  appropriate  turbulence  models 
are  seldom  incorporated  into  the  modeling  effort  mainly  because  of  the  additional 
difficulties  involved.  The  poor  understanding  of  turbulence  characteristics  in  break¬ 
ing  waves  seriously  impedes  the  choice  of  right  turbulence  model.  For  example, 
Miyata  (1986)  tried  to  simulate  wave  breaking  using  the  surface  segment  method 
without  including  any  turbulence  model.  Wang  and  Su  (1993)  also  attempted  to 
simulate  the  wave  breaking  on  sloping  beaches  using  the  VOF  method  without 
including  any  turbulence  model.  In  their  models,  the  energy  dissipation  due  to  the 
turbulence  is  represented  by  the  numerical  dissipation  partly  rather  than  the  ap¬ 
propriate  turbulence  model.  Such  simulations  do  not  provide  complete  information 
for  breaking  waves  because  the  important  physics  of  turbulence  is  not  explored  in 
the  computation. 

In  order  to  account  for  the  turbulence  effect  on  the  mean  flow,  two  major  tur- 


3 


bulence  closure  models  are  available  which  can  in  principle  describe  approximately 
the  turbulence  generation,  transport,  and  deca}^  The  most  complete  model  is  the 
Reynolds  stress  closure  model,  which  attempts  to  close  the  Reynolds  stresses  trans¬ 
port  equation  directly  {Launder  et  ah,  1975).  Due  to  the  lack  of  real  understanding 
of  turbulence  behavior,  the  closure  problem  for  high  order  correlation  is  always  a. 
challenging  subject.  This  is  especially  true  for  the  pressure-strain  rate  correlation 
term,  which  redistributes  the  turbulence  energy  among  different  directions.  In  fact, 
based  on  different  assumptions,  there  are  currently  at  least  five  different  closure 
models  for  the  pressure-strain  rate  correlation  term  and  three  different  models  for 
the  diffusion  term  (Demuren  and  Sarkar,  1993).  The  validation  and  verification  of 
these  models  for  different  turbulent  flows  itself  is  a  very  difficult  task. 

The  k  —  e  model  is  another  important  turbulence  closure  model  on  the  lower- 
level  than  the  Reynolds  stress  closure  model  (Launder  et  ah,  1972;  Launder  and 
Spalding,  1972;  Launder  and  Spalding,  1974;  Rodi,  1980).  In  this  model,  instead  of 
seeking  the  direct  closure  of  the  Reynolds  stresses  transport  equations,  the  eddy- 
viscosity  assumption  is  made,  which  relates  the  Reynolds  stresses  to  k,  turbulence 
kinetic  energy,  e,  the  dissipation  rate  of  k,  and  the  rates  of  strain  of  the  mean  flow. 
The  value  of  k  and  e  can  be  obtained  by  solving  the  transport  equations  derived 
from  the  NSE.  The  assumption  of  eddy-viscosity  not  only  reduces  the  number  of 
transport  equations  (six  for  the  Reynolds  stresses  closure  model  but  only  one  for 
the  k  equation)  but  also  simplifies  the  closure  assumption.  However,  we  must 
realize  that  such  an  assumption  also  limits  the  application  of  the  model  to  the 
nearly  isotropically  diffusive  turbulent  flow.  At  the  current  stage,  we  decide  to  use 
k  —  €  model  as  a  preliminary  tool.  It  is  our  plan  to  further  implement  Reynolds 
stress  closure  model  in  the  future  when  more  turbulence  characteristics  in  breaking 
wave  is  understood  through  the  experimental  and  numerical  studies. 

In  this  report  we  will  review  a  numerical  model  which  calculates  the  ensemble 
averaged  flow  field.  The  free  surface  is  tracked  with  the  use  of  the  VOF  method. 
The  effects  of  turbulence  are  included  by  coupling  the  mean  momentum  equa¬ 
tions  and  the  transport  equations  of  turbulence  kinetic  energy  and  the  dissipation 
rate.  This  approach  will  allow  one  to  examine  the  free  surface  configuration,  mean 
velocity  and  pressure,  and  the  turbulent  intensity  in  broken  waves. 


2  Mathematical  Formulation 

In  this  section,  the  mathematical  equations  governing  the  flow  motion  and  turbu¬ 
lence  transport  and  dissipation  will  be  summarized  respectively. 
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2.1  Flow  Motion  Formulation 

2.1.1  Navier-Stokes  Equations  and  Boundary  Conditions 

The  motions  of  an  incompressible  fluid  can  be  described  by  the  Navier-Stokes 
equations  in  a  bounded  domain  Q: 


dui 


+ 


dui 

dxi 

dui 

dxj 


pdxi^^'^pdx, 


(1) 

(2) 


where  ij  =  1,2,3  for  three-dimensional  flows.  The  Navier-Stokes  equations  rep¬ 
resent  the  conservation  of  mass  and  momentum  per  unit  mass  in  which  u,-  denotes 
the  ?’-th  component  of  the  velocity  vector,  p  the  density,  p  the  pressure,  gi  the  i-th 
component  of  the  gravitational  acceleration,  and  r^j  the  viscous  stress  tensor.  For 
a  Newtonian  fluid,  =  2paij  with  p  being  the  viscosity  and  <J,j  = 
the  rate  of  strain  tensor.  Two  types  of  boundary  are  considered:  rigid  boundary 
Tr  and  free-surface  boundary  Vj.  Along  the  rigid  boundary  the  velocity  of  the 
boundary,  Ui,  is  prescribed  and  the  no-slip  boundary  condition  requires 


ui  =  Ui  (3) 

On  the  free  surface  the  continuity  of  stress  components  must  be  required.  De¬ 
noting  n  as  the  unit  normal  on  the  free  surface  and  rn  as  the  projection  of  n  on 
the  Xi,  the  continuity  of  the  normal  stress  component  can  be  expressed  as 


p-p 


dui  duj 
dxj  ^  dxi 


TiiUj  =  r„ 


(4) 


in  which  is  a  prescribed  normal  stress  applied  on  the  free  surface  and  the  surface 
tension  has  been  ignored.  For  three-dimensional  problems,  two  unit  tangential 
vectors,  =  1,2),  are  needed  to  define  the  local  tangent  plane  on  the  free 

surface  and  is  defined  as  the  projection  of  on  the  x,.  If  external  tangential 
stress  components  are  applied  on  the  free  surface,  r/  and  the  continuity  of 
tangential  stress  components  becomes 


'  dui  duj 
^dxj  dxi 


flit  j 


fc  =  1,; 


(5) 


In  addition  to  the  stress  continuity  boundary  condition  on  the  free  surface, 
which  is  also  referred  to  as  the  dynamic  free  surface  boundary  condition,  the  kine¬ 
matic  boundary  condition  must  also  be  satisfied.  Before  a  wave  breaks,  such  a 


boundary  condition  assumes  the  velocity  continuity  that  ensures  the  free  surface 
to  be  material  surface.  Since  a  material  surface  always  consists  of  the  same  parti¬ 
cles,  the  total  derivative  of  any  physical  property  associated  with  the  free  surface 
particles,  which  is  expressed  as  G{x,t),  must  vanish  on  the  free  surface.  Thus, 


dt  dxi 


-0 


(6) 


After  a  wave  breaks,  the  free  surface  is  no  long  a  material  surface  and  equation 
(6)  is  not  valid.  In  the  present  model,  in  conjunction  with  the  volume  of  fluid 
method,  the  G  function  is  chosen  as  the  density  function,  i.e.,  G{x,t)  =  p{x,t). 
Equation  (6)  then  becomes  a  generally  true  statement  for  an  incompressible  fluid, 
before  and  after  the  wave  breaks.  If  the  density  remains  a  constant,  equation  (6) 
is  satisfied  automatically.  However,  in  the  volume  of  fluid  method,  the  density  of 
a  fluid  in  each  computational  cell  is  defined  as  the  averaged  density  in  the  cell. 
Therefore,  in  the  computational  cell  where  the  free  surface  appears,  the  “density” 
of  the  fluid  is  less  than  the  real  density  of  the  fluid,  while  the  “density”  of  the  fluid 
in  the  computational  cell  occupied  by  air  is  zero.  By  tracking  the  “density”  change 
in  each  computational  cell,  one  is  able  to  estimate  the  motion  of  the  free  surface. 


2.1.2  Reynolds  Equations  and  Boundary  Conditions 

Equations  (1)  -  (6)  describe  the  governing  equation  and  boundary  conditions  for  a 
wide  range  of  flow  motions  of  an  incompressible  fluid  with  a  free  surface,  including 
potential  flows,  laminar  and  turbulent  flows.  Once  the  suitable  initial  conditions 
for  the  velocity  and  pressure  field  and  the  boundary  configurations  are  specified, 
the  formulation  is  completed.  These  equations  can  be  solved  directly  by  numerical 
methods.  However,  in  the  case  of  turbulent  flows  with  high  Reynolds  number 
(Re),  the  resolution  for  small  scale  turbulent  fluctuations  is  so  high  that  the  direct 
numerical  simulation  (DNS)  is  extremely  difficult. 

An  alternative  way  to  understand  turbulent  flows  is  to  examine  the  mean  (en¬ 
semble  averaged)  flow  field  with  the  consideration  of  the  influence  of  turbulence. 
Both  the  velocity  field  and  the  pressure  field  can  be  split  into  mean  component 
and  turlsulent  fluctuations  as  follows: 


Ui  =  {Ui)  +  u'i 

p  =  (p)  +  p 

1  _  _1_ 

p  {p)  p 


(7) 

(8) 

(9) 


in  which  (  )  denotes  the  mean  quantities  and  the  prime  “'”  represents  the  turbulent 
fluctuations.  Therefore,  (u-)  =  (p)  =  {y)  =  0.  By  substituting  (7),  (8),  and  (9) 
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into  (1)  and  (2)  and  taking  the  ensemble  average  of  the  resulting  equations,  one 
obtains  the  governing  equations  for  the  mean  flow  field 


d{ui) 

dxi 


d{ui) 

dt 


djuj) 

dxj 


=  0 


d{p)  ,  ^  1 

"1  9i  I 


{p)  dx.  ■  ■  (p)  dxj 


-{ 


p  dxi  p  dxj 


djujUj) 

dxj 


(10) 


(11) 


These  equations  are  also  called  the  Reynolds  equations.  It  is  noted  that  com¬ 
pared  with  the  classical  Reynolds  equations,  there  are  two  additional  terms  con¬ 
tributed  by  the  correlations  between  the  density  fluctuation  and  the  gradient  of 
pressure  fluctuation  and  viscous  stress  fluctuation.  Since  we  assume  the  constant 
density  within  the  fluid,  these  two  terms  only  become  significant  near  the  free  sur¬ 
face.  At  this  time,  the  importance  of  these  correlations  is  totally  unclear  to  us  and 
thus  we  choose  to  temporarily  neglect  them  in  the  following  computations.  There¬ 
fore,  the  only  factor  that  will  be  taken  into  account  in  the  mean  flow  computation 
is  the  Reynolds  stress. 


R.,  =  -ipKu'u])  (12) 

Since  all  the  solid  boundary  conditions  and  dynamic  free  surface  boundary 
conditions  are  linear  in  terms  of  the  velocity  and  pressure,  equations  (3),  (4),  and 

(5)  apply  also  to  the  mean  flow  field  except  that  all  the  instantaneous  quantities 
are  replaced  by  the  corresponding  ensemble  averaged  quantities.  As  for  equation 

(6) ,  since  there  exist  nonlinear  terms  in  the  advection  of  function  G{x,t)  =  p, 
the  additional  terms  will  appear  in  the  mean  flow  equation.  After  performing  the 
ensemble  average  for  equation  (6)  using  the  similar  procedure  to  that  for  the  mass 
and  momentum  conservation  equations  except  for  letting  p  =  (p)  +  p\  we  obtain: 

The  right  hand  side  of  equation  (13)  represents  the  correlation  between  density 
fluctuations  near  the  free  surface  and  the  velocity  fluctuations.  Again,  because  of 
the  poor  understanding  of  such  correlation,  we  will  temporarily  neglect  it.  Equation 
(13)  becomes  the  same  as  equation  (6)  except  that  the  mean  quantities  replace  the 
instantaneous  quantities. 

It  is  noted  that  for  the  turbulent  flow,  the  free  surface  of  the  mean  flow  is  no 
longer  clearly  defined.  The  fluctuations  of  the  free  surface  caused  by  the  velocity 
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fluctuations  in  the  neighborhood  will  generate  a  region  with  the  finite  thickness 
and  the  variable  mean  density  {p)  that  ranges  from  pj  to  0  (air).  This  region  can 
be  regarded  as  the  mean  free  surface  region.  The  thickness  of  this  region  is  entirely 
dependent  upon  the  turbulence  intensity. 

Realizing  the  extreme  difficulties  of  the  closure  problem  for  the  fluctuations 
of  the  free  surface,  we  shall  neglect  this  effect  by  assuming  the  free  surface  of 
the  mean  flow  will  be  clearly  defined  even  in  the  turbulent  flow,  which  suggests 
(p)  =  p  =  pj  throughout  the  fluid.  Thus,  from  here  on,  the  (  )  sign  for  the 
mean  density  will  be  dropped.  The  above  approximation  enables  us  to  apply 
the  constant  density  assumption  within  the  well-defined  computational  domain  to 
simplify  the  formulation.  For  example,  the  viscous  term  in  equation  (11)  can  now 
be  rewritten  as  g^(2i/((Tjj))  with  u  —  mIp  being  the  kinematic  viscosity.  In  the 
following  derivation  of  A:  -  e  model,  it  is  this  modified  viscous  term  that  will  be 
used  rather  than  the  original  one  in  equation  (11). 

2.2  Turbulence  Transport  Model 

To  solve  Reynolds  equations  for  the  mean  flow,  one  must  relate  Reynolds  stresses 
to  the  mean  velocity.  Extensive  research  work  has  been  done  to  seek  the  proper 
closure  model  for  the  Reynolds  stresses  in  the  1970s  (e.g.  Launder  et  al.  1972, 
1975).  One  of  the  most  successful  approaches  is  the  fc  -  c  model  in  which  k  is  the 
turbulent  kinetic  energy  and  c  the  turbulent  dissipation  rate,  defined  as  follows 

k  =  ^{UiUi')  ,  c  =  ^  )  (14) 

The  basic  concept  of  the  k  —  t  closure  model  is  to  specify  the  relation  between 
Reynolds  stresses  and  the  rates  of  strain  of  the  mean  flow  as 

(«-u/)  =  -2ut{<Tij)  ^kSij  (15) 

in  which  Ut  is  the  eddy  viscosity,  depending  on  the  local  state  of  turbulence,  and 
Sij  is  the  Kronecker  delta.  The  eddy  viscosity  is  considered  as  a  function  of  k  and 
c.  Therefore,  the  remaining  tasks  are;  (1)  the  derivation  of  governing  equations 
for  k  and  e,  and  (2)  the  specification  of  the  relationship  among  lyf,  k  and  e. 

2.2.1  Transport  Equation  for  Turbulent  Kinetic  Energy,  k 

Multiplying  the  momentum  equations  for  the  total  velocity  (2)  by  the  turbulent 
velocity  component,  u  •,  and  taking  the  ensemble  average  of  the  resulting  equations, 
one  obtains  a  transport  equation  for  the  turbulent  kinetic  energy,  k: 
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dk  ,  ,  dk 


(16) 


The  left-hand  side  of  the  equation  calculates  the  rate  of  change  of  the  turbulent 
kinetic  energy  following  the  mean  flow  field.  The  first  term  on  the  right-hand  side 
represents  turbulent  diffusion  of  k  through  the  turbulent  pressure  work  done  and 
the  fluxes  of  kinetic  energy.  The  second  term  on  the  right-hand  side  represents  the 
molecular  diffusion  of  k,  which  is  usually  much  smaller  than  the  turbulent  diffusion. 
The  third  term  on  the  right-hand  side  denotes  the  rate  of  change  of  k  due  to  the 
working  of  the  Reynolds  stresses  against  the  mean  flow  gradients.  Therefore,  this 
term  represents  the  exchanges  between  the  mean  flow  energj^  and  the  turbulent 
flow  energy  and  it  is  also  called  the  production  term  P.  The  last  term  in  (16)  is 
the  energy  dissipation  rate  e  caused  by  the  viscous  stress. 

The  governing  equation  for  k  as  written  in  (16)  still  cannot  be  solved  because 
higher  order  correlations,  (u-p)  and  are  introduced.  In  principle  one  could 

derive  governing  equations  for  these  higher  order  correlations.  Unfortunately,  even 
higher  correlations  will  appear  in  the  resulting  equations.  Certain  closure  con¬ 
ditions  must  be  imposed  to  define  the  problem.  Following  Launder  k  Spalding 
(1972)  and  Rodi  (1980),  the  turbulent  diffusion  term  in  (16)  is  modeled  using  a 
gradient-diffusion  hypothesis: 


(17) 


in  which  cr^  is  an  empirical  diffusion  constant.  Substituting  (14),  (15)  and  (17)  in 
(16),  one  obtains  an  equation  for  turbulent  kinetic  energy,  k, 


dk  ,  y  dk 


dxj 


dxj 


( — I" 

Ok  OXj 


pl't 


d{ui) 


dxj  '  dxi  J  dxj 


—  e 


(18) 


in  which  the  continuity  equation  for  the  mean  flow,  (10),  has  been  employed. 
To  solve  the  above  equation,  information  on  the  mean  flow  field,  (u,),  which  is 
described  by  the  Reynolds  equations,  (11),  and  the  turbulent  dissipation  rate  e  is 
required.  Therefore,  one  needs  to  derive  another  equation  for  e. 


2.2.2  The  Transport  Equation  for  Turbulent  Dissipation  Rate,  e 

The  e  equation  can  be  derived  by  differentiating  the  Xi  component  of  the  momen¬ 
tum  equation  (2),  multiplying  the  resulting  equation  by  udujdxj,  and  taking  the 
ensemble  average  of  the  resulting  equation.  Thus, 
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^  dt  ^  du[du[du'  ^  /  dhi', 
dt^  dxj  dxkdxjdxj  \  dxkdxj^ 


dxi 


'  dxj  dxj 


dujdu- 

)+{  '  ' 


2  o  ^  ^  ^ 

I 

d{ui) 


dt 


pu-. 


dxi 


dxidxk  \  dxk 


,  du\  d'^{xH) 
M^k—) 


dxj  dxkdxj 


(19) 


Similar  to  the  exact  ^-equation,  (16),  the  exact  e-equation  contains  many  higher 
order  correlations.  The  physical  meaning  of  each  term  on  the  right-hand  side  of 
the  equation  can  be  given  as  follows:  The  first  term  represents  the  production 
by  vortex  stretching  due  to  the  turbulent  vorticity.  The  second  term  represents 
the  viscous  dissipation  due  to  the  spatial  gradients  of  turbulent  vorticity.  The 
third  term,  which  is  in  the  form  of  spatial  divergence,  represents  the  molecular  and 
turbulent  diffusion  of  e.  The  last  two  terms  on  the  right-hand  side  of  (19)  represent 
the  production  due  to  the  interaction  between  turbulent  correlations  and  the  mean 
velocity  gradients. 

To  close  the  problem  the  right-hand  side  terms  of  (19)  will  be  approximated. 
First  the  diffusion  terms  are  modeled  by  using  the  gradient-diffusion  assumption, 
i.e. 


•  r  d  (ut  dt\ 

Diffusion  oi  e  ^  -jr —  I  — — 
axj  \a^axjj 

where  is  an  empirical  constant.  Secondly,  the  difference  between  the  production 
and  dissipation  will  be  modeled  as  (Rodi  1980) 


(Production  -  Dissipation)^  — )■  Cu—i't 

in  which  Cit  ^■nd  C2t  are  two  more  empirical  coefficients  related  to  production  and 
dissipation  of  e,  respectively.  It  is  also  noted  that  the  kj e  denotes  the  characteristic 
decay  time.  Equation  (19)  can  now  be  simplified  to  be 


d{ui)  d{uj)\  d{ui)  _n  ^ 
dxj  dxi  )  dxj  k 


de  ,  \  de 
m  +  '  dx, 


d  ( ut  de 


dxi 


+  C\^—Ut 


'd{ui)  ^  d{uj) 
^  dxj  dxi 


d{u,) 

dxj 


C2.J  (20) 


From  the  dimensional  analysis,  the  eddy  viscosity  has  a  unit  of  (length^/time). 
Thus,  an  empirical  relation  among  r'f,  k  and  e  can  be  written  as 
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where  Cd  is  an  empirical  constant. 

Although  the  closure  assumptions  employed  in  the  k  —  t  equations  are  crude, 
this  model  has  been  successfully  used  to  predict  many  complex  flows.  The  empirical 
coefficients  appeared  in  the  model  are  surprisingly  universal.  The  recommended 
values  for  these  coefficients  are  (Rodi  1980): 


Cd  =  0.09,  Cl.  =  1.44,  Ca.  =  1.92,  =  1.0,  a.  =  1.3  (22) 

2.2.3  Determinations  of  the  coefficients  in  k  —  e  model 

Although  the  coefficients  for  the  k  -  e  model  are  somehow  regarded  to  Idc  stan¬ 
dard  for  many  types  of  turbulent  flows,  it  is  still  worthwhile  to  know  how  these 
coefficients  have  been  determined  and  what  assumptions  have  been  made.  These 
assumptions  also  suggest  the  limitations  of  the  model  and  the  possibilities  for  fur¬ 
ther  improvement. 

It  is  realized  that  the  original  k  —  e  equations  are  derived  from  the  basic  flow 
motion  equations,  i.e.,  Navier-Stokes  equations  and  Reynolds  stresses  equations, 
with  no  additional  assumptions.  The  k  —  e  model,  however,  consists  of  a  lot  as¬ 
sumptions  in  order  to  close  the  high  order  correlations.  Besides  those  assumptions 
in  the  turbulence  model,  we  also  assume  that  the  fundamental  turbulence  char¬ 
acteristics  does  not  change  under  different  flow  fields.  For  example,  we  assume 
that  the  same  k-e  model  with  the  same  coefficients  should  apply  equally  to  both 
well-developed  turbulent  flows  and  transient  turbulent  flows.  The  major  reason 
to  make  such  assumption  is  that  many  turbulence  measurements  that  are  used  to 
estimate  the  model  parameters  can  only  be  performed  and  analyzed  under  simple 
and  quasi-steady  state  (well-developed  turbulence).  While  most  engineering  ap¬ 
plications  of  turbulent  flows  are  complicated  and  transient.  The  gap  between  the 
theoretical  study  of  a  simple  turbulent  flow  and  the  practical  application  of  the 
complicated  turbulent  flows  has  not  been  resolved.  In  our  study,  we  retain  the 
above  assumption  to  study  the  complicated  and  transient  turbulent  flow  with  the 
use  of  turbulence  model  whose  coefficients  are  essentially  determined  under  the 
simple  situations. 

2.2. 3.1  Determination  of  €2^ 

To  determine  C2t,  we  will  consider  the  simple  dissipative  turbulent  flow.  In 
such  flow,  there  is  no  mean  velocity  gradient  and  thus  there  is  no  turbulence 
production.  The  typical  example  is  one-dimensional  grid  turbulence  in  x  direction 
that  is  created  by  making  flow  pass  certain  grid  system.  If  the  advection  is  much 


stronger  than  diffusion,  Chr  would  be  the  only  coefficient  appearing  in  equations 
(18)  and  (20).  For  the  well  developed  turbulence,  we  also  neglect  the  time  derivative 
term  so  that  k  and  e  are  functions  of  x  only.  The  simplified  equations  read: 


(«)^  =  -e  (23) 

clx 

(u)^  =  -cX  (24) 

ax  k 

The  general  solution  for  above  equations  system  are; 


k  —  ko{l  +  ^x)  ”  (25) 

c  =  eo(l+ea^)-”-‘  (26) 

with  ^  and  n  coefficients  and 

C.  =  (27) 

n 

The  measurement  of  the  decay  rate  of  k  in  the  grid  turbulence,  -n,  indicates  that 
C2c  lies  in  the  range  of  1.8  and  2.0  (Rodi,  1980).  In  most  models,  is  selected 
to  be  1.9. 

2. 2. 3.2  Determination  of  Cd 

Next,  we  turn  to  estimate  the  value  of  Cdi  the  coefficient  that  relates  the  value 
of  k  and  t  to  eddy  viscosity  Ut.  For  local  equilibrium  shear  layers  that  generally 
exist  in  constant-stress  turbulent  boundary  layers,  P  =  e.  Because  in  the  shear 
layer,  P  =  where  y  is  the  coordinate  normal  to  the  mean  flow 

(u),  and  according  to  the  definition,  e  =  Cd^i  we  have: 


C'rf  = 


(28) 


According  to  the  measurements,  the  value  of  ^  in  such  flows  is  around  0.3  that 
gives  Cd  —  0.09. 


2. 2. 3. 3  Determination  of  Cu 

In  the  log-law  region  of  a  fully  developed  turbulent  channel  flow,  it  is  assumed 
that  P  «  e  and  all  the  advection  is  negligible.  The  mean  velocity  gradient  can  be 
expressed  as  ^  ^  with  y  again  being  the  coordinate  normal  to  the  mean  flow 


direction,  being  the  frictional  velocity,  and  k.  =  0.41  being  the  von 

Karman  constant.  r„,  is  the  cros,s-stream  shear  stress  on  the  wall.  It  is  ready  to 
show  that  in  such  a  log-law  region,  the  k  remains  constant  over  y,  k  —  ^ 

and  III  =  —{u'v).  We  will  discuss  this  in  more  detail  later  when  we  propose  the 
wall  boundary  condition  for  the  k  —  t  model.  With  all  the  assumptions  mentioned 
above,  the  equation  for  e  (20)  can  be  reduced  to; 


Cjk'^  (Pt 
(j,e  dy"^ 


/ d^ ,2 


(29) 


By  substituting  the  solution  form  for  e  into  the  above  equation,  we  obtain: 


Cu  =  (30) 

As  long  as  the  coefficient  Cd,  C'2e,  and  are  known,  the  value  of  Cu  is  fixed.  At 
this  time,  we  have  already  determined  the  value  of  Cd  and  C'2e,  the  remaining  task 
is  then  to  determine  the  value  of  cr^  so  that  we  can  estimate  Cic- 

2.2. 3.4  Determinations  of  and  ak 

It  is  expected  that  the  value  of  cr^  and  ak  are  close  to  unity  that  implies  that 
both  k  and  e  diffuse  at  roughly  the  same  rate  as  the  mean  velocity.  By  the  computer 
optimization  (Launder  et  ah,  1972),  the  value  of  a^  and  ak  are  chosen  to  be  1.3 
and  1.0,  respectively.  The  value  of  (Te  fixes  the  value  of  Cu  to  be  1.49.  However, 
in  most  applications  of  the  k  —  e  model,  this  coefficient  is  adapted  as  1.44  as  first 
appearing  in  Launder  and  Spalding  (1974). 

2.2.4  Boundary  Conditions  for  k  and  e 

Although  the  boundary  conditions  for  the  total  flow  field  along  the  solid  boundaries 
and  the  free  surface  have  been  described  in  (3),  (4),  and  (5),  additional  boundary 
conditions  must  be  given  in  terms  of  k  and  e.  Since  the  physics  of  breaking  waves 
are  not  well  understood,  especially  near  the  free  surface,  it  is  extremely  difficult  to 
prescribe  boundary  conditions  for  k  and  e.  For  the  situations  where  the  external 
forces  are  absent  on  the  free  surface  it  seems  to  be  reasonable  to  assume  that 
turbulence  does  not  diffuse  across  the  free  surface.  Consequently,  the  normal  flux 
of  k  and  e  should  vanish  on  the  free  surface. 


dxi 


Hi  —  0 


de 

dxi 


Tli 


(31) 


In  theory,  the  turbulence  vanishes  at  wall  which  dictates  that  both  k  and  c 
become  zero  at  wall.  However,  in  most  of  practical  computations,  the  resolution 
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cannot  be  so  small  to  resolve  the  viscous  sublayer.  Thus,  the  boundary  conditions 
for  k  and  e  are  generally  specified  in  the  turbulent  boundary  layer  instead  of  right  on 
the  wall.  In  the  turbulent  boundary  layer  the  cross-stream  shear  stress  dominates 
and  remains  a  constant.  Invoking  the  boundary  layer  approximation,  we  have; 


{uv) 

dy  dy'^ 


=  0 


(32) 


where  y  is  the  coordinate  normal  to  the  wall.  By  taking  the  integration  from  the 
wall  to  the  place  out  of  the  viscous  layer  where  the  viscous  effect  can  be  neglected, 
we  have: 


d{u) 


I  yzzO 


=z  == 


(33) 


dy  p 

As  discussed  in  section  2.2.3. 3,  the  mean  velocity  gradient  in  this  region  can  be 
expressed  as: 


d{u)  _ 
dy  Ky 


(34) 


Integrating  the  above  equation,  we  have  so  called  logarithmic-law  profile  for  stream- 
wise  velocity: 


M  =  ife  (e^)  (35) 

n*  K  \  u  J 

where  E  =  9.0  for  a  smooth  wall.  Because  the  dissipation  rate  is  approximately 
the  same  as  production  rate,  e  =  P,  from  (33)  and  (34)  we  have: 


n; 


/  / 


From  (15)  the  eddy  viscosity  I't  can  be  obtained 


(36) 


(itv) 

dy 


(37) 


Namely,  the  eddy  viscosity  is  proportional  to  the  distance  from  the  wall  in  the 
turbulent  boundary  layer.  Substituting  (36)  and  (37)  into  (21)  yields 


Equations  (36)  and  (38)  constitute  the  boundarj'  conditions  for  k  and  e  at  the 
computational  points  immediately  adjacent  to  the  solid  boundary.  The  frictional 
velocity  can  be  found  from  (35)  once  the  velocity  field  has  been  calculated. 


2.3  Summary  of  Governing  Equations 

Before  we  leave  this  section,  it  is  useful  to  have  a  brief  summary  of  the  assumptions 
made  so  far  in  deriving  the  Reynolds  equations  for  mean  flow  motions  and  the  k  —  t 
model  for  turbulence  transport.  The  final  formula  that  will  be  used  in  the  later 
numerical  model  will  also  be  presented.  Since  only  the  mean  quantities  are  involved 
in  the  Reynolds  equations  and  k-t  model,  the  symbol  for  the  ensemble  average, 
(  ),  will  be  dropped  from  herein  for  simplicity.  By  invoking  the  eddy  viscosity 
assumption  (15)  and  neglecting  density  fluctuations  near  free  surface,  equations 
(10)  and  (11)  that  governing  the  mass  and  momentum  conservation  of  mean  flow 
can  be  rewritten  as: 


dui 

dxi 

dui  dui 
dxj 


=  0 


1  dp 


dr, 


dt 


p  dxi  Qx 


(39) 

(40) 


The  density  change  equation  that  can  be  used  to  track  the  free  surface  is  written 


as: 


dp  ^  a 


(41) 


The  k  —  e  model  reads: 


dk  dk 
dt 


d 


dxi  dxj 


(^t  .  dk 

CTk  dx 


3i 


dt 

dt  dxj 


dxj  [ 


.  .dt 


dui 


dui  „  t‘ 


-a 


(42) 

(43) 


and  the  total  stresses  r;,  in  equation  (40)  is  related  to  k,  c,  and  aij  in  the  following 


way. 


P 

2(z/  +  Cd~)<yij 


Equations  (39)  to  (43)  will  serve  as  the  governing  equations  to  be  solved  numeri¬ 
cally. 
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Numerical  Model 


In  this  section,  we  will  review  the  numerical  methods  used  to  solve  the  turbu¬ 
lent  flow  problem  as  described  in  the  previous  section.  Currently,  we  will  restrict 
ourselves  to  the  two-dimensional  problems  only. 


3.1  Schematics  of  Computational  Domain 

The  finite  difference  method  will  be  used  throughout  the  computation.  At  the 
beginning,  the  whole  computational  domain  is  discretized  by  the  mxn  rectangular 
cells  as  sketched  in  figure  1.  All  scalar  quantities,  i.e.,  p,  k,  e,  the  volume  of  fluid 
(VOF)  function  F,  and  the  openness  function  9c,  are  defined  in  the  center  of  the 
cells.  The  first  three  scalars  have  been  introduced  before,  and  the  last  two  scalars 
will  be  defined  later.  The  vector  and  vector-related  quantities,  i.e.,  the  x-  and 
y-components  of  the  mean  velocities,  u,  v,  and  the  openness  functions  on  the  cell 
faces.  Or  and  Ot,  are  defined  in  the  cell  faces  as  shown  in  figure  1.  Again,  the  last 
two  functions  will  be  defined  in  the  following  text. 

As  mentioned  before,  the  exact  location  of  the  free  surface  will  not  be  pursuited 
in  the  present  model.  Instead,  the  density  change  in  each  cell  will  be  tracked  so 
that  the  location  of  the  free  surface  can  be  identified.  We  assume  that  the  density 
is  a  constant  pj  for  the  fluid,  and  the  averaged  density  in  the  cell  is  defined  as, 
p  =  4^,  where  Vf  is  the  volume  of  the  fluid  in  the  cell  and  K  is  the  volume  of  the 
air  in  the  cell.  Before  solving  the  density  change  equation  (41),  we  further  simplify 
the  problem  by  using  the  normalized  averaged  density  or  the  VOF  function,  F  = 
pjpf.  With  the  aid  of  the  VOF  function,  we  may  identify  the  cell  types  readily. 
For  example,  the  empty  (E)  cell  is  defined  as  the  cell  with  F  =  0,  where  the 
computation  will  be  skipped;  the  surface  (S)  cell  is  defined  as  the  cell  with  F  >  0 
and  adjacent  to  at  least  one  empty  cell,  where  the  free  surface  boundary  condition 
will  be  applied;  and  the  interior  (I)  cell  is  defined  as  the  cell  which  has  no  empty 
neighbour  cells,  where  the  main  computation  will  be  conducted. 

Another  important  issue  is  how  to  represent  the  arbitrary  shape  of  solid  bound¬ 
ary.  In  other  models,  the  solid  boundary  is  generally  redefined  as  sawtooth-shape 
surface  to  fit  the  cell  boundary  (Lemos,  1992).  Based  on  our  numerical  tests,  such 
a  treatment  may  generate  the  significant  spurious  reflection  during  the  wave  shoal¬ 
ing.  In  the  current  model,  the  more  flexible  way  to  represent  the  solid  boundary, 
partial  cell  treatment,  is  adapted.  In  such  treatment,  the  solid  object  is  modeled 
as  the  special  case  of  the  two-phase  flow  with  infinite  density.  The  openness  func¬ 
tions  are  introduced  at  the  cell  itself  and  on  the  cell  faces.  At  the  cell  center,  Oc 
is  defined  as  the  ratio  of  space  not  occupied  by  the  solid  object  (thus  open  to  the 
fluid)  to  the  whole  cell  area.  On  the  cell  faces.  Or  {0i)  is  defined  as  the  length  open 
to  the  fluid  to  the  whole  length  of  the  right  (top)  cell  boundary.  Therefore,  similar 
to  the  VOF  function,  the  openness  functions  can  provide  the  information  whether 


the  cell  is  the  solid  object  or  obstacle  (0),  the  fluid  (air)-solid  boundary  (FA-0), 
or  the  fluid  (air)  domain  (FA).  The  only  difference  between  the  VOF  function'and 
the  openness  functions  is  that  the  former  is  a  time- varying  function  and  the  latter 
are  not. 

Figure  1  illustrates  the  definition  of  different  cell  types  based  on  the  information 
of  the  VOF  function  and  openness  functions.  In  the  following  conte.xt,  the  imple¬ 
mentation  of  the  VOF  function  advection  and  the  treatment  of  solid  boundary 
with  openness  functions  will  be  discussed  in  more  detail. 

3.2  Two-step  Projection  Method 

The  two-step  projection  method  (Chorin  1968,  1969)  has  been  used  to  solve  the 
Reynolds  equations.  The  first  step  is  to  introduce  an  intermediate  velocity  u,-, 
which  carries  the  correct  vorticity,  as 


-  < 


At 


=  —u 


^du'j 

n _ J_ 

dxi 


+  di  T 


dxj 


(44) 


in  which  the  superscript  indicates  the  time  level  and  At  is  the  time  step  size.  The 
subscripts  =  1,2  for  two-dimensional  problems,  i.e.,  Ui  —  u,  U2  =  u,  r’l  =  a;, 
and  yi  =  y.  Therefore  (44)  is  the  forward  time  difference  equation  of  the  Reynolds 
equation  without  the  pressure  gradient  term.  The  intermediate  velocity  u,  does 
not,  in  general,  satisfy  the  continuity  equation.  If  all  information  on  the  ?r-th  time 
level  is  known,  (44)  is  an  explicit  formula  for 

The  second  step  is  to  project  the  intermediate  velocity  field  onto  a  divergence- 
free  plane  to  obtain  the  final  velocity: 


At 

dxi 


p"  dxi 


=  0 


(45) 

(46) 


One  can  combine  (44)  and  (45)  to  show  that  the  Reynolds  equations  are  satisfied 
approximately,  with  pressure  gradient  being  evaluated  at  the  {n  -f-  l)-th  time  level. 
The  continuity  equation  is  satisfied  by  (46).  Taking  the  divergence  of  (45)  and 
applying  (46)  to  the  resulting  equation  yields 
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which  is  also  called  the  Pois.son  Pressure  Equation  (PPE). 


3.3  Spatial  Discretization  in  Finite  Difference  Form 

In  the  two-step  projection  method  the  spatial  derivations  of  the  velocih^  compo¬ 
nents  and  the  pressure  field  need  to  be  expressed  in  finite-difference  forms.  As 
discussed  before  and  shown  in  figure  1,  the  present  scheme  calculates  the  velocity 
components,  u  and  v,  on  the  vertical  and  horizontal  cell  faces,  respectively.  The 
pressure  and  other  scalars  such  as  k,  e,  and  the  volume  of  fluid  function  F  are  de¬ 
fined  at  the  cell  center.  It  is  noted  that  in  the  finite  difference  form,  some  variables 
a.re  needed  at  the  place  where  they  are  not  originally  defined.  For  instance,  u  at 
the  horizontal  cell  face  and  v  at  the  vertical  cell  face.  In  such  circumstances,  the 
linear  interpolation  will  be  used.  The  most  commonly  used  interpolated  variables 
are  given  as  follows: 
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As  appearing  in  equation  (44),  all  the  advection  terms  and  diffusion  terms 
will  be  evaluated  at  the  n-th  time  step.  The  advection  terms  in  the  x-momentum 
equation  +  are  evaluated  at  the  front  (right)  face  of  the  cell.  The  advection 

terms  in  the  y-momentum  equation  are  calculated  at  the  top  face  of  the 

cell  (Figure  2).  Thus 
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In  equation  (54)  one  needs  to  specify  the  spatial  derivatives  of  the  horizontal  veloc¬ 
ity.  To  calculate  these  spatial  derivations  of  the  horizontal  velocity  component  in 
(54),  one  can  use  the  combination  of  the  upwind  scheme  and  the  central  difference 
scheme.  For  instance 
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where 

Axa  =  Axi+i  +  Axi  +  a  sgn  {ui+i,j)  (Ax^+i  -  Axi)  (60) 

Similarly 
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In  the  above  equations,  the  coefficient  a  is  the  weighting  factor  between  the 
upwind  method  and  the  central  difference  method.  When  a  =  0,  the  finite  differ¬ 
ence  form  becomes  the  central  difference,  while  a  =  1,  the  finite  difference  form 
becomes  the  upwind  difference.  In  practice,  a  is  generally  selected  in  the  range  of 
0.3  to  0.5  to  produce  the  stable  and  accurate  results.  The  finite  difference  form  for 
the  advection  terms  in  the  y-momentum  equation  (55)  can  be  similailj'^  obtained. 

After  neglecting  the  term  -jkSij  in  the  total  stresses,  which  can  be  lumped  into 
the  pressure  term,  the  gradients  of  the  total  stresses  in  (44)  can  be  represented  by 
the  diffusion  terms  that  are  written  as 
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for  the  x-momentum  equation  and 
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for  the  y-momentum  equation.  Once  again  the  diffusion  of  the  momentum  in 
the  x-direction  is  calculated  at  the  right  face  of  the  computational  cell,  while  the 
diffusion  in  the  y-direction  is  computed  at  the  top  face  of  the  cell.  The  central 
differencing  method  is  used  to  express  the  derivatives.  Thus,  the  first  term  of  (65) 
can  be  written  in  the  following  finite  differences  from 
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The  second  term  in  (65)  can  be  written  as 
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in  which 
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Similar  finite-difFerence  formulas  can  be  obtained  for  the  diffusion  terms  in  the 
j/-momentum  equation  (66). 

In  the  second  step  of  the  projection  method,  the  Poisson  Pressure  equation  (47) 
needs  to  be  solved.  In  the  two-dimensional  form  (47)  can  be  rewritten  as 
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Since  the  pressure  is  evaluated  at  the  center  of  the  computational  cell,  the  Poisson 
equation  is  also  evaluated  at  the  center.  Thus,  the  left-hand  side  of  the  equation 
can  be  expressed  as: 
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The  right-hand  side  of  (69)  is  expressed  as 
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which  are  the  known  quantities.  In  (70a)  and  (70b)  the  density  on  the  side  of  the 
cells  can  be  obtained  by  the  linear  interpolation,  i.e. 


Pi+k,j  ~ 


P”, Ax,-+i  -h  pr+i,j Ax,- 
Ax,-  -f  Ax,-+i 


(71) 


Substituting  (70)  and  (71)  into  the  Poisson  equation  yields  a  set  of  linear  alge¬ 
braic  equations  for  the  pressure  field  that  can  be  solved  by  standard  matrix  solvers. 
In  our  model,  the  conjugate  gradient  method  with  the  preconditioner  of  incomplete 
Cholesky  decomposition  is  used  to  solve  the  resulting  sparse  and  symmetric  system 
of  equations. 
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3.4  Volume  of  Fluid  Method 

As  explained  at  the  beginning  of  this  section,  the  volume  of  fluid  method  is  a 
means  to  identify  different  types  of  computational  cells  and  thus  can  be  used  to 
track  the  free  surface  motion.  Letting  p{x,  y,  t)  =  F{x,  y,  t)pj  and  substituting  this 
definition  and  equation  (39)  into  equation  (41),  we  obtain  the  transport  equation 
for  F(.T,y,f), 


For  a  computational  cell  centered  at  (f,  j)  the  above  equation  can  be  rewritten  in 
the  following  finite-diflference  form 
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in  which  Fg,  F£,  Ff  and  F^  denote  the  F  values  on  the  right,  left,  top  and  bottom 
faces  of  the  computational  cell,  respectively.  Assigning  the  proper  F  values  on 
each  face  of  the  cell  requires  some  information  of  the  free  surface  configuration. 
For  instance,  the  F  values  in  the  computational  cells  centered  at  {i,j)  as  shown  in 
figure  3a  and  36  are  the  same.  However,  the  free  surface  configurations  are  quite 
different.  Even  if  the  velocities  on  the  right  face  of  the  computational  cells  are  the 
same,  one  would  expect  that  more  fluids  are  convected  across  the  face  in  case  (b) 
than  in  case  (a)  for  small  At.  For  this  reason  a  simple  free  surface  reconstruction 
scheme  is  needed.  In  this  report,  we  summarize  the  Hirt-Nichols  algorithm  in  which 
the  free  surface  is  reconstructed  either  horizontally  or  vertically  in  each  free-surface 
cell  based  on  the  F  values  at  the  n-th  time  step.  For  the  cases  shown  in  figure 
3,  the  reconstructed  free  surface  profiles  are  quite  different.  The  general  rule  is  as 
follows 
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Figure  3:  Two  types  of  free  surface  reconstruction  in  VOF 

The  first  derivatives  of  F  can  be  evaluated  by  using  the  central  difference 
scheme.  For  instance 
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in  which  is  the  average  value  of  F  at  the  cell  (i  +  l,i)  using  three  vertical 

neighboring  cells,  i.e., 
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Similar  formulas  can  be  obtained  for  and  F^-\^j- 

Once  the  free  surface  orientation  is  determined,  the  donor-acceptor  method  is 
used  to  advect  the  F  function.  For  simplicity  only  the  flux  in  the  x-dii'ection  will 


be  discussed.  Moreover,  the  .^-component  of  the  new  velocity  on  the  right  face 
of  the  computational  cell  is  assumed  to  be  positive,  i.e.,  ^  >  0.  In  the 

donor-acceptor  method,  the  cell  {i,j)  is  called  the  donor  cell,  while  the  cell  (f  +  l,j) 
is  called  the  acceptor  cell.  If  .  <  0,  the  role  of  the  donor  and  acceptor  cell  is 

switched. 

When  the  free  surface  configuration  in  the  donor  cell  is  determined  to  be  a 
horizontal  surface,  the  value  for  F  on  the  right  face  is  assigned  to  be  the  same  as 
that  in  the  donor  cell,  i.e. 
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On  the  other  hand  if  the  free  surface  configuration  in  the  donor  cell  is  a  vertical 
surface,  the  Fr  value  should  be  more  strongly  influenced  by  the  F  value  in  the 
acceptor  cell,  i.e. 
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However,  this  acceptor  method  may  create  the  so-called  over-filling  problem  in 
some  situations.  Consider  the  simple  one-dimensional  case  as  shown  in  case  (a) 
of  figure  3.  The  F  values  in  three  adjacent  cells  are  jP/Iij  =  1-0,  F-^j  =  0.8  and 
F.”  =0.  In  other  words,  there  is  a  sharp  front  in  the  {i,j)  cell.  If  the  acceptor 

method  is  used,  equation  (77)  indicates  that  Fr  =  0. 

From  the  one-dimensional  version  of  (73),  one  finds  that 
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in  which  F"  =  0.8  and  FV  —  1.0.  Therefore,  if  the  Courant  number  Cr  = 
u^^ljAtfAx  is  greater  than  0.2,  the  F  value  at  the  (n  -f  1)  time  step  will  be 
greater  than  one.  To  overcome  this  “over-filling”  problem,  a  corrector  term  is 
introduced  in  the  acceptor  method,  i.e. 
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in  which  denotes  the  value  of  F  at  the  cell  upstream  of  the  donor  cell,  which 
is  equal  to  id—\  when  >  0.  K  is  a  factor  used  to  prevent  spurious  advection 

of  a  small  amount  of  fluid  and  is  often  set  to  be  0.1. 


3.5  The  k  —  e  equations 

The  same  nonuniform  finite-difference  mesh  is  used  to  discretize  the  k—t  equations, 
(42)  and  (43).  Both  k  and  e  are  defined  at  the  center  of  a  computational  cell. 
Symbolically,  the  k  —  e  equations  in  the  finite-difference  form  can  be  written  as 
(Lemos,  1992); 
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In  (81)  and  (82)  the  advection  terms  {FkX,  FkY,  FeX  and  FeY),  viscous  diffusion 
terms  {VISk  and  VISe),  and  the  production  terms  Pij  need  to  be  defined  and 
described  in  space.  Since  the  expressions  of  these  terms  for  the  k  equation  are  very 
similar  to  those  for  the  t  equation,  only  the  details  for  the  k  equation  are  given 
here. 

First,  we  examine  the  advection  terms.  The  x-component  of  the  advection  term 
can  be  written  as 
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Note  that  the  derivative  of  k  is  evaluated  at  the  n-th  time  level,  which  makes  the 
scheme  explicit.  The  finite  difference  form  of  the  spatial  derivative  for  k,  dkfdx, 
is  the  same  as  that  for  dufdx.  Thus 
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where 


MV  ^KizBzhL 

Aii_i 

M'j  ^  -  kr.i 

V9^/,+i,,  M+i 


(85) 

(86) 


and  7  is  the  coefficient  weighting  the  relative  importance  between  the  central  differ- 
®  ence  and  the  upwind  scheme,  similar  to  a  in  the  finite  difference  form  of  momentum 

advection.  In  the  actual  computation,  7  is  generally  chosen  to  be  the  unity  that 
makes  the  advection  terms  discretized  by  the  upwind  scheme.  Such  scheme  ensures 
the  stable  solution  for  the  k  equation  under  any  circumstance.  However,  it  is  found 
that  for  most  computations,  as  long  as  7  is  greater  than  0.5,  almost  the  same  stable 

•  solutions  will  be  obtained. 

The  y-component  of  the  advection  term  can  be  expressed  in  the  similar  form 

as  (83) 
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in  which 
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The  viscous  term  in  (81)  is  defined  as: 
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The  first  term  in  above  equation  can  be  discretized  with  the  central  difference 
method, 


/A 

dk 

X  - 

Kf'k  ^  A'+lJ  ' 

r*  f^Y 
1-Tj 

dx 

i.9  ^ 

Axi 

Axi 


-  Ki  C'-  I  .. 


i-hjJ 


Ax 


(^k 


Ax-  I 
*  2 


(92) 


and  the  second  term  can  be  similarly  obtained, 
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where  the  value  of  (^  +  ^)"  defined  at  the  cell  face  can  be  obtained  with  linear 
interpolation,  i.e., 


(-+■')”  = 
\Gk  J 
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The  production  term  P,j  is  evaluated  at  both  the  n-th  and  (n  +  l)-th  time 
steps.  At  the  n-th  time  step,  we  have; 


with  the  following  finite  difference  forms: 
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The  value  of  can  be  similarly  obtained. 
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3.6  Stability  Analysis  for  Numerical  Model 

The  governing  equations  (40)  to  (43)  are  essentially  the  transient  advection-diffusion 
equations.  The  explicit  finite  difference  form  for  these  equations  are  subject  to  the 
numerical  instability  unless  certain  stability  criteria  are  satisfied.  Due  to  the  intrin¬ 
sic  nonlinearity  in  the  governing  equations,  the  general  stability  analysis  method, 
i.e.,  von  Neumann  method,  is  not  applicable.  In  order  to  overcome  this  difficulty, 
we  first  linearize  the  governing  equations  and  then  perform  the  standard  von  Neu¬ 
mann  stability  analysis  to  obtain  the  stability  criteria  for  the  linear  case: 


At  <  min 


At  <  min 
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Ax  ’  Ay j 
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(100) 

(101) 


The  above  criteria  are  obtained  by  assuming  the  advection  term  is  discretized 
by  upwind  scheme.  The  first  criterion  arises  from  the  stability  requirement  for  the 
advection  term  and  the  second  from  the  diffusion  term.  In  practice,  due  to  the 
effect  of  nonlinearity  and  the  combination  of  the  central  difference  method  that 
will  further  restrict  the  stability  condition,  some  empirical  coefficients  that  are 
less  than  unity  are  multiplied  to  the  original  time  step  restriction  to  ensure  the 
practical  stability  in  the  computation.  These  coefficients  are  A  for  (100)  and  |  for 
(101). 


3.7  Computational  Cycle 

Finally,  the  complete  cycle  for  updating  the  field  variables  for  one  time  step  is 
summarized  as  follows: 

1.  Compute  the  tentative  velocities  using  equation  (44). 

2.  Apply  the  tentative  boundary  conditions  (tangential  stress  and  divergence  free) 

on  the  free  surface  as  well  as  the  inflow  boundary  condition  if  necessary. 
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3.  Update  the  pressure  field  according  to  equation  (47),  which  automatically 

incorporates  the  normal  stress  boundary  condition. 

4.  Obtain  the  final  velocities  according  to  equation  (45). 

5.  Apply  again  the  boundary  conditions  on  the  free  surface. 

6.  Update  k  a.nd  c  value  using  new  velocities. 

7.  Update  the  VOF  function. 

8.  Apply  the  final  boundary  condition,  mainly  to  the  newly  filled,  cell  that  was 

originally  empty  in  the  previous  time  step 

4  Detailed  Numerical  Treatments  in  RIPPLE  and 
Modifications  in  Present  Model 

In  this  section,  some  detailed  numerical  treatments  in  RIPPLE  will  be  discussed. 
A  few  modifications  have  been  made  to  improve  the  accuracy  of  the  numerical 
computations.  The  major  modifications  appear  on  the  applications  of  boundary 
conditions  on  the  free  surface.  Efforts  have  also  been  made  to  improve  the  accuracy 
of  numerical  solution  to  the  PPE  and  the  momentum  conservation  equation  near 
the  solid  boundaries. 


4.1  Boundary  Conditions  on  Free  Surface 

In  RIPPLE,  the  velocity  field  at  the  new  time  step  for  both  the  interior  cells  and 
surface  cells  are  obtained  by  using  the  projection  method.  Then  the  divergence  free 
criterion  is  imposed  on  the  free  surface  cell  to  adjust  the  velocity  at  the  interface 
between  the  surface  cell  and  empty  cell.  For  example,  for  cell  (i,  j),  Wj+ij, 
and  V-  ■  I  are  obtained  by  the  projection  method,  but  is  obtained  by  using  the 

divergence  free  criterion,  +  ||  =  0?  which  leads  to  ==  ~  “ 

Ui_y)  (Figure  4).  _  .  .  j  , 

In  RIPPLE,  there  is  no  explicit  dynamic  boundary  conditions  applied  on  the 
free  surface.  The  only  criterion  enforced  on  the  free  surface  is  divergence  free  that 
ensures  the  mass  conservation  on  the  free  surface.  It  is  found  that  this  approach 
will  create  large  spurious  velocity  on  the  free  surface  that  significantly  affects  the 
accuracy  of  velocity  and  vorticity  computation  nearby.  The  modifications  are  made 
in  the  following  aspects; 
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Figure  4:  Boundary  conditions  and  new-filled  cell  treatment  on  free  surface 


4.1.1  Tangential  Stress  Boundary  Conditions 


In  our  computations,  the  air  effect  is  neglected  so  that  the  tangential  stress  bound¬ 
ary  condition  becomes:  |^  + =  0-  For  nearly  horizontal  free  surface,  the  above 
equation  becomes:  =  0.  For  the  cell  {i  -  l,j),  conventionally  the  boundary 

condition  is  applied  at  point  {i  -  |,  j  +  |)  (Lemos,  1992).  Thus,  the  corresponding 


ihiii  =  0  (Figure  4).  The 


finite  difference  form  becomes:  - . 

outside  fictional  tangential  velocity  ni_i/2j+i  can  be  found  from  the  above  expres¬ 
sion  provided  that  the  other  three  quantities  are  known.  We  have  found  that  this 
treatment  can  still  generate  large  errors  in  the  velocity  computations  on  the  free 
surface.  The  errors  in  the  velocity  computation  are  mainly  due  to  the  fact  that  the 


pressure  computed  on  the  free  surface  usually  contains  large  errors.  The  errors  in 
the  pressure  computation  on  free  surface  cell  arise  from  the  fact  that  in  RIPPLE 
the  normal  stress  cannot  be  applied  onto  the  exact  free  surface  in  general  when 
the  PPE  is  solved.  Fortunately,  the  pressure  computed  just  one  cell  away  from  the 
free  surface  is  very  accurate.  This  phenomenon  is  caused  by  the  particular  method 
in  RIPPLE  to  solve  the  PPE  and  will  be  discussed  later  in  detail. 

Realizing  that  the  pressure  computed  on  the  free  surface  may  contain  large  er¬ 
rors,  which  in  turn  affects  the  tangential  velocity  computations,  we  decide  to  apply 
the  tangential  stress  Ijoundary  condition  one  cell  down  from  the  free  surface.  For 
instance,  for  cell  (i  -  l,i),  instead  of  applying  the  boundary  condition  at  the  point 
“  2 ’ 2)’  ^PP^y  boundary  condition  at  the  point  (f  —  i  —  |)  a-nd  obtain 
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by  using  and  By  so  doing,  we  obtain  the  velocities 

between  surface  cells  not  by  using  the  projection  method,  but  by  employing  the 
tangential  stress  boundary  condition.  Because  the  computed  pressure  of  the  cell 
beneath  the  free  surface  is  very  accurate,  the  new  way  for  obtaining  the  velocity 
on  the  free  surface  is  much  more  accurate  than  the  conventional  method. 

4.1.2  Newly  Filled  Cell  Treatment 

After  updating  the  VOF  function,  the  new  surface  cells  may  be  created,  which 
were  originally  empty  at  the  previous  time  step.  RIPPLE  fails  to  apply  the  correct 
boundary  condition  for  this  situation.  Chen  et  al.  (1995)  addressed  this  problem 
and  proposed  the  momentum-capturing  method  to  define  the  velocity  for  those 
newly  filled  cells.  The  basic  idea  for  their  method  is  to  assume  the  newly  filled  cell 
carries  the  momentum  of  its  neighborhood.  For  example,  in  figure  4,  cell  (f — 1,  j-fl) 
is  empty  at  the  end  of  n-th  time  step  but  is  partially  filled  at  the  end  of  {n  -1-  l)-th 
time  step.  In  RIPPLE,  =  0  at  the  end  of  (n  -f  l)-th  time  step.  This  will 

affect  the  accuracy  of  velocity  computations  near  the  free  surface  at  (n-|-2)-th  time 
step.  On  the  other  hand,  Chen  et  al.  (1995)  assigned  ^be  end 

of  {n  4-  l)-th  time  step  that  ensures  the  accuracy  of  the  velocity  computation  on 
the  free  surface.  Our  model  follows  the  same  idea  as  Chen  et  al.’s  to  define  the 
velocity  at  the  newly  filled  cells. 

4.2  PPE  Computation 

During  the  computation  of  PPE,  the  normal  stress  boundary  condition  is  naturally 
incorporated.  In  most  of  situations,  the  air  effect  and  the  normal  stress  due  to  the 
fluid  on  the  free  surface  is  neglected,  which  simplifies  the  normal  stress  condition 
to  be  p  =  0  on  the  free  surface  from  equation  (4).  Since  the  free  surface  seldom 
passes  the  cell  center  where  the  pressure  is  defined,  the  application  of  the  normal 
stress  boundary  condition  on  the  free  surface  becomes  a  challenging  task.  In  the 
MAC  method,  the  accuracy  of  applying  normal  stress  boundary  condition  can  be 
improved  by  using  the  irregular  star  method  (Chan  and  Street,  1970)  and  the  micro 
cells  treatment  (Raad  et  al.,  1994).  These  methods,  however,  cannot  be  readily 
extended  to  the  VOF  method  where  the  free  surface  location  is  not  exactly  defined. 

Realizing  this  limitation,  RIPPLE  uses  another  approach  to  treat  the  normal 
stress  boundary  condition.  Instead  of  attempting  to  apply  the  pressure  on  the  exact 
location  of  the  free  surface,  the  approximation  is  made  to  assume  the  free  surface 
cell  is  the  cell  with  variable  density  fluid  in  which  the  density  at  the  cell  center  is 
represented  by  the  VOF  function,  pc  =  pjF  (Figure  5).  Such  an  approximation 
is  consistent  with  the  approximation  we  made  before  in  equations  (6),  (13),  and 
(41)  that  the  free  surface  has  a  finite  thickness  filled  by  stratified  fluids.  Thus, 
the  free  surface  is  always  assumed  to  be  located  at  the  center  of  empty  cell  just 


above  the  free  surface  cell  during  the  PPE  computation.  The  above  approximation, 
though  crude,  is  simple  to  use  and  provides  the  very  accurate  pressure  solution  one 
cell  awa)^  from  the  free  surface  cell.  The  disadvantage  of  the  method  is  that  the 
pressure  computed  right  in  the  surface  cell  may  contain  large  errors.  These  errois 
may  contaminate  the  tangential  velocities  computation  on  the  free  surface. 

Computed  p  using 
projection  method 
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Figure  5:  An  example  of  PPE  computations  for  one-dimensional  problem 


To  understand  why  the  pressure  computed  on  the  free  surface  may  contain  large 
errors  and  the  pressure  one  cell  below  is  almost  exactly  computed,  we  will  examine 
the  following  simple  example.  The  one-dimensional  case  with  the  still  water  free 
surface  above  the  fourth  cell  center  (Figure  5)  will  be  examined.  The  constant  grid 
size  Ay  will  be  used  in  the  numerical  computation.  Assuming  initially  there  is  no 
flow  motion,  we  expect  that  the  fluid  will  remain  still  and  the  pressuie  distribution 
should  be  hydrostatic.  Numerically,  we  hope  that  by  advancing  one  time  step 
from  initially  still  condition,  the  fluid  can  stay  still  and  the  pressure  hydrostatic. 
Following  the  standard  project  method  described  before,  we  will  first  find  out 
that  the  tentative  velocity  without  using  the  pressure  information.  Then  we  solve 
the  PPE  with  known  tentatively  velocity  divergence.  By  imposing  the  boundary 
condition  on  the  free  surface  ps  =  0  and  on  the  solid  bottom  pi  —  p2-,  it  is  very 
easy  to  find  the  solution  to  the  PPE  to  be:  p4  =  p^igAy,  p3  =  [p^i-t- p^i)gAy,  and 
P2  =  (P41  +  pgi  -t- P2i)gAy.  Now,  let’s  inspect  the  pressure  in  the  cell  right  beneath 
the  free  surface,  cell  3.  Since  p  is  represented  by  the  VOF  function  p  =  Fpj  and 
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the  value  of  F  at  the  cell  face  can  be  obtained  by  linear  interpolation  of  F  in  the 
cell  center,  we  may  rewrite  =  pjg[{F4  +  |F3)Aj/].  It  is  ready  to  show  that  the 
pressure  computed  there  is  exactly  hydrostatic  and  so  are  all  cells  beneath.  The 
pressure  on  the  free  surface  is  found  to  be  p4  =  ^F4pjgAy.  It  is  ready  to  observe 
that  the  pressure  computed  on  the  free  surface  cell  may  deviate  from  the  true 
pressure  and  thus  can  affect  the  final  tangential  velocity  computation  for  the  two- 
dimensional  problem.  One  effective  way  to  reduce  the  influence  of  pressure  errors 
is  to  adapt  the  new  method  to  apply  the  tangential  stress  boundary  condition  on 
the  free  surface  as  introduced  in  the  previous  section. 

During  the  testing  of  RIPPLE,  another  error  is  found  in  RIPPLE  which  reduces 
the  accuracy  of  solution  to  the  PPE  when  a  variable  mesh  size  is  used.  Referring 
back  to  equation  (70a),  when  the  spatial  derivative  of  p  at  the  cell  face,  for  example, 

is  discretized,  we  will  have  .  Instead  of  using  + 

Arci+i),  which  represents  a  first  order  accurate  discretization  for  the  variable  grid, 
RIPPLE  uses  =  |(^  +  a^)’  zeroth  order  accurate. 

*+ 2 

4.2.1  Partial  Cell  Treatment 

Partial  cell  treatment  is  a  special  tool  used  in  RIPPLE  to  handle  the  unstructured 
interior  obstacle  and  solid  boundary.  As  introduced  at  the  beginning  of  previous 
section,  the  basic  idea  behind  this  technique  is  that  the  obstacle  can  be  mod¬ 
eled  as  a  special  case  of  two-phase  flow  with  infinite  density.  Compared  with  the 
conventional  way  to  treat  irregular  solid  wall  that  creates  the  saw-tooth  shape  of 
boundary,  partial  cell  treatment  partially  block  the  cell  face  and  cell  itself  according 
to  the  real  geometry  of  the  boundary.  With  the  use  of  the  partial  cell  treatment, 
the  variables  in  cells  or  cell  faces  are  redefined  as  the  multiplication  of  openness 
coefficients  with  original  variables,  representing  the  mean  values  there.  Near  the 
obstacle,  the  openness  coefficients  are  less  than  unit  that  makes  the  mean  quan¬ 
tities  smaller  than  their  original  values.  The  governing  equations  for  the  general 
turbulent  flow  can  be  rewritten  as: 


d(6ui) 

dt 


d{0ui) 
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-j-  Ouj 


d{0Ui) 
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3  dxi 


(102) 

(103) 


In  RIPPLE,  the  advection  terms  is  written  as  for  the  conservative 

scheme  with  one  9  missing.  It  is  found  that  RIPPLE’s  approach  sometimes  causes 
instability  during  the  runup  when  6  is  small  in  some  location  where  Ui  is  relatively 
large,  while  equation  (103)  eliminates  such  numerical  instability. 
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5  Numerical  Examples 

In  this  section,  we  will  verify  our  model  by  performing  the  numerical  experiments 
for  a  series  of  testing  problems.  These  testing  problems  are  in  the  order  of  in¬ 
creasing  complexity.  The  first  group  of  the  numerical  experiments  is  for  laminar- 
flow  without  turbulence.  The  purpose  of  such  testing  is  to  demonstrate  the  capa¬ 
bility  of  the  solver  to  the  NSE.  Without  losing  generality,  we  use  three  examples 
for  numerical  testing.  One  is  the  motion  of  initially  tilting  fluid  in  a  rectangular 
tank.  The  numerical  results  are  compared  to  the  linear  wave  solution.  The  second 
example  is  the  solitary  wave  propagation  on  the  constant  water  depth.  The  nu¬ 
merical  results  of  the  wave  profiles  are  compared  to  the  analytical  solutions  based 
on  the  Boussinesq  equation.  The  third  example  is  the  nonbreaking  solitary  wave 
runup  on  a  uniform  slope.  The  numerical  results  are  compared  to  the  BIEM  results 
based  on  potential  flow  theory  and  experimental  measurements.  This  test  shows 
the  capability  of  the  partial  cell  treatment  in  handling  the  unstructured  bottom 
topography.  The  second  group  of  the  numerical  experiments  are  for  turbulence 
flows.  The  purpose  of  such  testing  is  to  validate  the  turbulence  model  used  in  the 
computation.  The  validation  of  the  turbulence  model  is  generally  very  difficult. 
For  a  well-developed  turbulent  flow  without  free  surface,  the  k  —  e  model  has  been 
sufficiently  verified  in  the  previous  studies  as  summarized  by  Rodi  (1980).  The 
verification  of  the  k  —  t  model  for  a  transient  flow  with  free  surface,  however,  is 
rarely  conducted  before.  The  example  we  choose  is  the  breaking  solitary  wave 
runup  on  beach  (Synolakis,  1986;  Zelt,  1991).  The  simulated  wave  profiles  will  be 
compared  with  experimental  measurements. 


5.1  Laminar  Oscillating  Fluids  in  a  Tank 

In  this  example,  flow  oscillations  in  a  rectangular  tank  are  simulated.  The  free 
surface  of  the  fluid  in  the  tank  is  initially  at  an  inclined  position  (Figure  6(a)). 
After  the  initial  moment,  the  fluid  begins  to  move  under  gravity.  Since  the  flow 
motions  are  restricted  in  the  tank  on  both  sides,  the  motions  can  be  decomposed 
into  infinite  modes  of  standing  waves  with  the  maximum  wave  length  of  L  =  2W, 
where  W  is  the  width  of  the  tank.  Mathematically,  the  free  surface  displacement 
r]{x,t)  can  be  written  as: 


T]{x,  t)  =  ^  o„  sin(A;„x)  cos(a;„t)  (104) 

n=l 

where  A'„  =  ^  is  the  wave  number  of  the  7r-th  mode,  the  frequency  of  the  n-th 

mode  and  can  be  determined  by  the  dispersion  relationship,  u>n  =  \/gkn  tanh(A:„d), 
with  d  being  the  still  water  depth.  a„,  the  wave  amplitude  of  the  n-th  mode,  can 
be  determined  by  using  the  initial  condition. 
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(105) 


_  sL 

”  2n^7r^ 

where  is  the  slope  of  the  initially  inclined  free  surface.  The  analytical  solution 
is  accurate  as  long  as  the  linear  wave  theory  is  valid.  In  the  actual  computa¬ 
tions  of  the  analytical  solution,  40  modes  are  used  which  are  sufficient  to  resolve 
the  the  free  surface  variations.  In  the  numerical  computations,  the  width  of  the 
tank  IT  =  0.277r,  the  water  depth  d  =  0.04m,  and  the  initial  free  surface  slope 
s  =  0.1.  Total  40  X  80  meshes  (non-uniform  in  vertical  direction)  are  placed  onto 
the  computational  domain.  The  time  step  At  is  adjusted  automatically  during 
the  computation  so  that  Cr  =  max{^,  ^)  =  0.3.  This  guarantees  the  practi¬ 
cal  stability  while  at  the  same  time  maintains  the  computational  efficiency.  The 
molecular  viscosity  is  neglected  in  the  computation  to  simulate  the  potential  flow. 
It  has  been  shown  that  such  simplification  has  a  negligible  effect  to  the  numerical 
results. 

Figures  6(a)  to  6(1)  give  the  comparison  of  the  analytical  solution  (104)  and 
(105)  to  the  numerical  results.  The  comparisons  show  that  the  numerical  results 
agree  very  well  with  the  analytical  solution  after  about  one  period  for  the  slowest 
mode,  Ti  =  —  =  0.507s.  The  small  discrepancies  in  the  comparisons  could  be 
caused  by  both  the  linear  wave  approximation  for  finite  amplitude  waves  and  the 
numerical  truncation  errors  in  the  finite  difference  scheme.  The  VOF  method 
behaves  very  well  in  tracking  free  surface,  even  though  there  are  more  than  one 
major  modes  during  the  motion.  Such  characteristics  of  numerical  scheme  is  very 
useful  for  the  breaking  wave  simulation,  where  shorter  waves  are  generated  during 
the  breaking  process. 

It  is  noted  that  different  grid  systems  varying  from  20  x  30  to  80  x  120  have 
been  used  to  check  the  sensitivity  of  numerical  results  to  grid  system  used.  It  is 
found  that  the  numerical  results  computed  with  different  grid  system  are  fairly 
close.  In  the  next  examples,  the  similar  sensitivity  analysis  for  grid  systems  has 
also  been  performed  and  the  final  choice  of  mesh  size  is  determined  by  considering 
both  computational  efficiency  and  numerical  accuracy. 

5.2  Solitary  Waves  Propagation  on  Constant  Water  Depth 

The  solitary  wave  propagation  on  the  constant  water  depth  is  a  classical  testing 
problem  for  a  wave  simulation  model.  The  solitary  wave  is  a  finite  amplitude  wave 
with  the  permanent  shape.  The  nonlinearity  and  frequency  dispersion  are  perfectly 
balanced  during  the  wave  propagation.  The  analytical  solution  for  the  wave  profile 
can  be  derived  from  the  Boussinesq  equation  (Lee  et  ah,  1982), 


nn  „  .  .  , 

4sm(-^)  —  2sm(777r) 
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where  a  is  the  wave  amplitude  and  c  =  ^gd(l  +  j)  is  the  wave  celerity. 

Three  cases  are  investigated  in  this  example,  with  different  ratios  of  the  wave 
amplitude  to  the  still  water  depth  of  a/d  =  0.1,  ajd  =  0.2,  and  ajd  =  0.3  The 
still  water  depth  is  d  =  Im  for  three  tests.  In  the  numerical  computations,  the 
same  discretizations  of  Ax  =  O.bn  and  Ay  =  0.0l77r  are  used  for  three  tests 
with  0?7z  <  X  <  100m.  The  time  step  At  is  automatically  adjusted  during  the 
computation  to  satisfy  the  stability  constraints  by  both  the  advection  and  diffusion 
processes,  i.e.,  equations  (100)  and  (101).  The  velocity  u  and  v  and  free  surface 
displacement  y  are  specified  at  the  left  boundary  based  on  the  Boussinesq  analytical 
solution  for  a  solitary  wave  (Lee  at  al,  1982). 

Figure  7  shows  the  comparisons  between  the  simulated  wave  profiles  and  the 
analytical  solutions.  In  Figure  7(a),  the  case  of  ajd  =  0.1  is  presented.  It  is 
observed  that  the  numerical  results  agree  with  the  analytical  solution  perfectly  even 
after  the  propagation  of  about  lOOd.  This  test  verifies  the  accuracy  of  the  numerical 
model  under  small  value  of  a/d.  Figures  7(b)  and  7(c)  shows  the  case  of  ajd  =  0.2 
and  a/d  =  0.3.  It  is  found  that  the  discrepancies  between  the  numerical  results 
and  the  analytical  solutions  increase  as  the  increase  of  a/d,  especially  after  a  long 
distance  of  propagation.  These  discrepancies  could  be  caused  by  two  factors.  One  is 
that  as  the  increase  of  a/d,  the  numerical  dissipation  increases  accordingly  because 
the  numerical  dissipation  rate  is  proportional  to  (9(aAx),  while  u  is  proportional 
to  a/d.  On  the  other  hand,  as  the  increase  of  a/d,  the  Boussinesq  equation,  which 
neglects  the  higher  order  nonlinear  and  dispersive  terms,  may  not  be  sufficient  to 
approximate  the  original  NSE. 

5.3  Nonbreaking  Solitary  Wave  Runup  on  Beach 

In  this  example,  the  nonbreaking  solitary  wave  runups  on  a  steep  beach  with  the 
slope  of  30°  are  investigated.  The  toe  of  the  beach  is  6.49m  away  from  the  wave 
maker.  The  still  water  depth  is  0.16m  and  the  wave  amplitude  is  0.027m.  Numer¬ 
ical  computations  using  BIEM  method  are  performed  for  the  whole  domain,  from 
the  wave-maker  to  the  shoreline.  The  number  of  nodes  is  121  on  the  free  surface, 
51  along  the  bottom,  11  on  the  sloping  beach,  and  6  on  the  wave  maker.  The 
time  step  is  At  =  0.01s.  The  numerical  computations  using  the  VOF  method  is 
performed  in  a  smaller  domain  near  the  shoreline  with  4.49m  <  x  <  6.99m  and 
-0.16?77  <  y  <  0.11m.  The  computational  domain  is  discretized  with  a  250  x  90 
uniform  grid  systems  with  Ax  =  0.01m  and  Ay  =  0.003m.  The  fixed  time  step 
At  =  0.004s  is  used  that  produces  a  stable  solution.  The  velocity  u  and  v  and 
free  surface  displacement  rj  are  specified  at  the  left  boundary  similarly  to  that  in 
section  5.2.  The  velocity  is  measured  by  a  particle  image  velocimetry  (PI\  )  system 
(Adrian,  1991).  Figure  8  gives  the  sketch  of  the  computational  domain  and  the 
measurement  section. 

The  results  of  the  nonbreaking  solitary  wave  are  shown  in  figure  9  which  includes 
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Figure  8:  Computational  domain  for  VOF  and  measurement  section  for  PIV 

the  VOF,  BIEM,  and  PIV  results  for  the  free  surface  profile  and  internal  velocity. 
Figure  9(a)  to  8f  show  the  time  series  of  solitary  wave  runup,  rundown,  and  re¬ 
runup  on  the  steep  beach.  The  overall  agreement  of  the  free  surface  profiles  between 
the  experimental  measurements  and  numerical  computations  are  excellent.  The 
comparisons  of  the  internal  velocity  at  different  time  are  also  very  good.  Some 
discrepancies  are  observed  between  experimental  data  and  numerical  results  when 
wave  runup  reaches  its  highest  point  and  rundown  reaches  its  lowest  point.  It 
is  noted  that  the  BIEM  results  and  the  VOF  method  results  always  agree  very 
well.  Since  the  BIEM  results  for  potential  flow  are  very  accurate,  they  can  serve 
here  as  the  reference  for  both  the  experimental  measurements  by  the  PIV  and  the 
numerical  results  by  the  VOF  method. 

Figure  9(a)  shows  the  snapshot  when  the  wave  is  climbing  up  the  slope  {t  = 
6.38s).  It  is  interesting  to  observe  that  the  particles  in  the  wave  front  move  in  the 
same  direction  as  the  beach  slope  orientation.  The  particles  at  the  farthest  runup 
tip  possess  the  maximum  velocity.  Velocity  decreases  downward  and  backward 
due  to  the  negative  pressure  by  the  runup  wave.  In  this  frame,  experimental  data 
and  numerical  results  agree  very  well.  Figure  9(b)  (<  =  6.58s)  shows  the  velocity 
distribution  when  the  wave  almost  reaches  its  highest  runup  point.  Except  for  the 
very  small  tip  that  still  has  the  forward  and  upward  velocity,  the  major  portion  of 
fluid  on  the  slope  has  started  the  rundown  process.  In  this  frame,  the  comparison 
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Figure  9:  Nonbreaking  solitary  wave  runup  on  a  sloping  beach 
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Figure  9:  Nonbreaking  solitary  wave  runup  on  a  sloping  beach 
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of  velocity  field  agrees  well  only  in  the  rear  part  of  wave.  Near  the  wave  front,  the 
discrepancies  between  experimental  data  and  numerical  results  are  obvious,  mainly 
due  to  different  velocity  directions.  Such  discrepancies  imply  the  existence  of  phase 
errors.  It  is  found  that  the  PIV  system  has  a  timing  error  in  the  range  of  0.01s. 
Such  a  small  timing  error  is  insignificant  normally,  considering  the  wave  period 
is  in  the  order  of  magnitude  of  Is.  However,  during  the  transition  from  runup 
to  rundown,  such  an  error  might  becomes  crucial  because  the  local  time  scale  for 
such  transition  process  is  in  the  order  of  0.1s  from  the  experimental  measurements. 
Thus,  we  expect  the  velocity  disagreement  in  the  wave  front  is  most  likely  caused 
by  the  timing  inaccuracy. 

Figure  9(c)  (t  =  6.78s)  shows  that  the  wave  is  running  down  on  the  beach  and 
figure  9(d)  {t  =  7.18s)  shows  that  the  wave  is  approaching  its  maximum  running 
down  point.  Both  the  wave  profile  and  velocity  field  comparisons  are  satisfactory 
in  these  two  frames.  After  wave  reaches  its  maximum  rundown  point,  the  positive 
pressure  starts  to  push  it  back  and  forms  the  secondary  runup.  Such  process  is 
shown  in  figure  9(d)  {t  =  7.38s)  and  figure  9(e)  {t  =  7.58s).  Small  discrepancies 
again  show  up  in  velocity  field  comparison  at  t  =  7.38s  when  the  wave  begins 
to  climb  up  the  slope.  Such  discrepancies,  which  is  expected  to  be  caused  by 
timing  inaccuracy  as  well  as  the  sensitivity  of  wave  during  transition  process,  is 
diminishing  as  time  precedes.  This  can  be  observed  in  figure  9(f). 


5.4  Breaking  Solitary  Waves  Runup  on  Beach 

In  this  example,  the  breaking  solitary  runup  on  the  beach  is  simulated.  The  detailed 
experimental  setup  for  this  testing  problem  is  referred  to  Synolakis  (1986)  and  Zelt 
(1991).  The  important  parameters  are  summarized  here.  The  still  water  depth 
varies  from  0.21m  to  0.29m.  The  beach  slope  is  2.88°  and  is  about  1/20.  A  solitary 
wave  which  has  a  wave  height  to  still  water  depth  ratio  of  0.28  was  generated. 
Wave  gauges  were  used  to  record  the  free  surface  displacement  during  the  runup 
and  rundown.  In  the  VOF  computation,  the  computational  domain  starts  from 
the  toe  of  beach  and  it  is  discretized  with  a  260  x  65  uniform  grid  system  with 
Ax  =  0.025m  and  Ay  =  0.005m.  The  time  step  is  automatically  adjusted  during 
the  computation  to  satisfy  the  stability  constraints  by  both  advection  and  diffusion 
processes,  i.e.,  equations  (100)  and  (101).  Similar  to  the  previous  nonbreaking 
solitary  wave,  both  the  velocities  and  free  surface  displacement  are  specified  on 
the  right  boundary  based  on  the  analytical  solution  for  solitary  wave.  Since  the 
k  —  c  model  is  employed  during  the  computations,  the  specification  of  the  inflow 
boundary  conditions  for  both  k  and  e  is  necessary.  As  shown  in  equation  (43), 
both  production  term  and  dissipation  term  for  e  become  singular  when  k  goes  to 
zero.  Thus,  in  general  low  levels  of  k  are  enforced  in  the  inflow  boundary,  i.e., 
k  =  with  Ut  =  SC  and  C  the  phase  velocity  of  wave  in  the  constant  water 
region.  In  our  computation,  S  is  chosen  to  be  around  0.003.  It  is  noted,  however. 


the  simulation  result  in  the  interior  region  is  basically  independent  of  the  value  of 
S  as  long  as  it  is  small  enough.  The  value  of  e  is  calculated  using  equation  (21),  for 
example,  e  =  Cd^  with  Vt  =  (v.  In  our  computation,  (  is  chosen  to  be  10.  Again, 
the  final  result  is  insensitive  to  the  choice  of  C  as  long  as  the  calculated  e  on  the 
inflow  boundary  remains  a  low  level. 

Since  the  water  depth  varies  in  the  experiment,  the  results  can  only  be  compared 
after  normalization.  Following  Zelt  (1991),  the  length  scale  is  normalized  by  the 
water  depth  and  the  time  scale  is  normalized  by  yg/d.  Figures  10(a)  to  10(h) 
show  the  comparisons  of  experimental  data  and  numerical  results  in  terms  of  the 
free  surface  profile.  The  whole  process  of  wave  breaking,  runup,  and  rundown  is 
provided.  The  agreement  between  the  simulated  wave  profiles  and  experimental 
data  is  excellent. 

Compared  with  the  numerical  results  of  Zelt  (1991)  obtained  from  the  Boussi- 
nesq  equation  and  with  those  of  Titov  and  Synolakis  (1995)  obtained  from  shallow 
water  equation,  our  model  provides  much  better  prediction,  especially  during  and 
after  wave  breaks.  Their  depth-averaged  models  tend  to  predict  larger  reflection 
which  is  not  observed  in  the  experiment  and  our  numerical  results  (Figure  10(d)). 
The  success  of  the  current  model  is  the  result  of  two  aspects.  First  is  that  our 
model  solves  the  Reynolds  equations  directly  and  thus  is  free  of  the  assumption  for 
pressure,  which  could  lead  to  significant  errors  during  and  after  the  wave  breaks. 
Another  is  that  the  more  elaborated  turbulence  model,  the  k  —  t  model,  provides 
more  realistic  turbulence  simulation  which  improves  the  free  surface  simulation. 

In  figure  11(a)  to  11(h),  the  magnitudes  of  the  turbulent  intensity  are  plotted. 
It  is  observed  that,  before  the  wave  breaks,  the  turbulence  is  mainly  concentrated 
in  a  very  thin  layer  near  the  bottom  and  the  quantities  of  those  turbulence  is 
relatively  small.  At  tJgJd  =  15,  wave  starts  to  break  and  the  increasing  large 
turbulence  is  generated  on  the  front  face  of  wave.  The  time  of  breaking  predicted 
by  the  numerical  model  is  consistent  with  what  suggested  by  Zelt  (1991).  In  the 
experiments,  the  breaking  process  occurs  between  tyjgjd  =  15  and  t-^fgjd  =  20. 
The  turbulence  generated  near  the  free  surface  is  left  over  as  the  broken  wave 
moves  towards  to  the  beach.  The  residual  turbulence  quickly  transport  to  the 
interior  region  where  it  starts  to  dissipate  because  of  the  small  mean  flow  velocity 
gradient.  This  can  be  observed  from  t^fgjd  -  20  and  =  25.  As  the  wave 

continues  to  climb  to  its  highest  point,  the  overall  turbulence  almost  disappears 
(ty^/d  =  45).  The  wave  then  starts  to  rundown  and  the  turbulence  energy  begins 
to  increase  again  due  to  the  backward  flow  motion.  The  intensity  of  turbulence, 
however,  is  relatively  smaller  than  that  generated  by  the  violent  breaking  process 
during  the  shoaling.  This  is  observed  at  t^jgjd  —  50. 
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Figure  11:  Simulated  turbulence  intensity  for  breaking  solitary  wave 
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Conclusions 


In  this  report,  we  have  given  a  detailed  review  of  the  numerical  model  simulating 
the  breaking  wave.  The  numerical  accuracy  and  the  capability  of  the  model  is 
verified  by  comparing  the  simulated  results  with  the  available  analytical  solutions 
and  experimental  data  for  both  laminar  flow  and  turbulent  flow.  It  is  found  that  for 
the  laminar  flow,  the  model  can  provide  very  accurate  results.  For  turbulent  flow 
such  as  breaking  wave,  the  comparison  for  the  wave  profile  is  excellent.  However, 
the  more  rigorous  validation  in  terms  of  turbulence  quantities  are  still  needed  for 
a  complete  verification. 
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